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Adequate simulation of non-adiabatic dynamics through conical intersection requires account for a non-trivial 
geometric phase (GP) emerging in electronic and nuclear wave-functions in the adiabatic representation. 
Popular mixed quantum-classical (MQC) methods, surface hopping and Ehrenfest, do not carry a nuclear 
wave-function to be able to incorporate the GP into nuclear dynamics. Surprisingly, the MQC methods 
reproduce ultra-fast interstate crossing dynamics generated with the exact quantum propagation so well as if 
they contained information about the GP. Using two-dimensional linear vibronic coupling models we unravel 
how the MQC methods can effectively mimic the most significant dynamical GP effects: 1) compensation for 
repulsive diagonal second order non-adiabatic couplings and 2) transfer enhancement for a fully cylindrically 


symmetric component of a nuclear distribution. 

I. INTRODUCTION 

The Born-Oppenheimer representation of the electron- 
nuclear wave-function introduces natural separation be¬ 
tween time/energy scales of electrons and nuclei in molec¬ 
ular systems. This separation allows one to consider nu¬ 
clear dynamics independently from that of the electronic 
subsystem reducing the number of involved degrees of 
freedom (DOF). This representation is uniquely defined 
through the electronic eigenvalue problem with fixed nu¬ 
clei and is conveniently available in numerous electronic 
structure packages. However, there are also a few com¬ 
plications associated with the inherent non-separability 
of dynamics in a general interacting many-body system. 
For example in many photochemical processes nuclear 
molecular dynamics cannot be adequately modelled on a 
single potential energy surface (PES) because for some 
nuclear configurations the separation between electronic 
PESs becomes comparable to the nuclear energy scale or 
even vanishes. The latter case often presents itself in the 
form of conical intersections (CIs). 1 ’ 2 Non-adiabatic dy¬ 
namics associated with such crossings not only results in 
transferring system population between electronic states 
but also in geometric phase (GP) effects. 3-8 The latter 
is caused by a nontrivial nuclear dependent geometric or 
Berry phase appearing in both nuclear Xj (R, t) and elec¬ 
tronic <f>j(r;R.) wave-functions within the adiabatic rep¬ 
resentation of the total electron-nuclear wave-function 9 

^(r,R,Q =X $ j(r; R )Xj(R,i)- (1) 

j 

Berry 10 and Longuet-Higgins 11 have shown that para¬ 
metric evolution of the electronic parts (r; R) around 
the point of eigenvalue degeneracies (the Cl seam) 
must change their signs, which makes <I>j(r;R) double¬ 
valued functions of nuclear DOF R. To preserve 
the single-valued character of the total wave-function, 
T(r, R, t), the nuclear part, Xj(R, t), must also have a 



FIG. 1. Low energy (a) and excited state (b) nuclear dynam¬ 
ics near CL 

double-valued character compensating that in the elec¬ 
tronic components. GPs in electronic wave-functions 
$7'(r; R) are needed to obtain non-adiabatic couplings 
(NACs) [($ !; (R)|V R ^(R))r and ($ 8 (R)| V^(R)> r ] 
correctly, 12 but NACs alone are not sufficient for simulat¬ 
ing nuclear dynamics properly. For the correct dynamics 
near CIs, the nuclear part must also have the GP result¬ 
ing in double-valued nuclear wave-functions y ? (R). Ig¬ 
noring the GP of nuclear wave-functions can lead to qual¬ 
itative distortion of non-adiabatic dynamics even in the 
the absence of a significant population transfer between 
crossing electronic states. 3-5,13 Interestingly, dynamical 
features associated with the GP are very different for low 
energy dynamics (Fig. la) and excited state dynamics 
(Fig. lb). The main manifestation of GP effects in low 
energy nuclear dynamics is destructive interference be¬ 
tween two paths around the Cl seam (Fig. la), 3,4,13 while 
in the excited state dynamics (Fig. lb), it is enhance¬ 
ment of population transfer for a fully cylindrical compo¬ 
nent of a nuclear wave-packet and compensation of a re¬ 
pulsive diagonal Born-Oppenheimer correction (DBOC), 
($*(R)| Vj l < f> i (R)) r . 6 The DBOC is usually neglected in 
non-adiabatic simulations for molecular systems based 
on its small value near the minimum of the ground state. 
However, at the intersecting manifold this term diverges 
to infinity (see Fig. 2), and its a priori neglect is not 
justified. 

One of the most popular ways to simulate non- 
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FIG. 2. Near CIs, the DBOC becomes very large and alters 
electronic PESs. In the absence of the GP, the DBOC can 
inhibit the access of a nuclear-wave packet to the CL 


adiabatic dynamics near CIs in large systems is us¬ 
ing mixed quantum-classical (MQC) approaches: surface 
hopping (SH) and Ehrenfest (EF) methods. 14 ’ 15 Unlike 
full quantum approaches, MQC methods propagate nu¬ 
clear DOF classically. As a result, they exhibit some 
well-known limitations of classical mechanics such as in¬ 
ability to model nuclear tunnelling and quantum inter¬ 
ference effects. The GP induced destructive interference 
in low energy dynamics (Fig. la) is a typical example of 
the latter effect. MQC methods have only the electronic 
part of the wave-function and thus cannot fully account 
for the GP, because even though the electronic function 
acquires the GP through parametric dependence on nu¬ 
clear evolution, nuclei evolve according to classical New¬ 
ton equations that do not have any GP contributions. 16 
In this context, it is quite surprising that short-time 
excited state deactivation dynamics (Fig. lb) of MQC 
methods agrees extremely well with that of the exact 
quantum propagation 1 '~ 19 for systems where GP effects 
were proven to be very influential. 6 

In this work we will explain how MQC methods em¬ 
ulate GP effects in Cl problems. We will restrict our 
attention to GP effects in excited state deactivation pro¬ 
cess (Fig. lb). As for low energy dynamics (Fig. la), SH 
and EF methods are not much better than a simple classi¬ 
cal dynamics because non-adiabatic transitions are well 
suppressed by the energy difference in areas accessible 
for classical trajectories. Currently, only the quantum- 
classical Liouvillc formalism has proven to be capable to 
capture GP effects in low energy dynamics. 20 

The rest of the paper is organized as follows. First, we 
introduce a diabatic two-dimensional linear vibronic cou¬ 
pling (LVC) model, although very simple representation 
of the Cl topology, this model has all quantities involved 
in quantum and MQC simulations in the analytical form. 
Thus, it allows us to compare various quantum and MQC 


methods on the same footing and to reveal the key com¬ 
ponents of the MQC schemes that are responsible for 
mimicking quantum GP effects. Second, to confirm our 
analysis we simulate non-adiabatic dynamics for a few 
molecular systems that provide a variety of dynamical 
regimes. Finally, we conclude the paper with a summary 
and further ramifications of our work. 


II. THEORY 
A. 2D LVC model 

The GP appears only in the adiabatic representation, 
however, it is more convenient to start with a model in 
the diabatic representation because a diabatic model will 
allow us to obtain the GP explicitly. 21 Moreover, the di¬ 
abatic representation can be exactly transformed to the 
adiabatic representation while the exact reverse transfor¬ 
mation does not exist in general. 22 Note that although 
diabatic wave-functions do not have GPs, simulations in 
the diabatic representation incorporate all GP effects im¬ 
plicitly and are considered to be exact. Thus, we start 
with introducing a diabatic model Hamiltonian 

ftD = r a „i 2+ (££), (2) 

where T 2 D = —h 2 2~ 1 (d 2 /dx 2 + d 2 /dy 2 ) is the nuclear ki¬ 
netic energy operator, x and y are mass-weighted nuclear 
coordinates, V\\ and V 22 are the diabatic potentials rep¬ 
resented by identical two-dimensional parabolas shifted 
in space and coupled by the U 12 potential 

V 11 = ^ [uJi(x + x 0 ) 2 + u%y 2 ] , V 12 = cy, (3) 

V 22 = ^ [wi(a; - x 0 f + uj%y 2 ] . (4) 

Here, uii s are harmonic frequencies for nuclear coordi¬ 
nates x and y, ilo are the minima of V\\ and V 22 po¬ 
tentials, and c is a coupling constant. Electronic DOF 
in H 2 t> are vectors |0i ia ) and |^>f la ) in a two-dimensional 
linear space. To obtain the corresponding adiabatic rep¬ 
resentation for the Hamiltonian H 2 y> one needs to diago¬ 
nalize the two-state potential matrix in Eq. (2) by a uni¬ 
tary transformation that defines the adiabatic electronic 
states as 

l^? di ) 
l^) 

where 

, 2 V 12 

0 = arctan —-—— (7) 

U 11 — V 22 

is a mixing angle between the diabatic states |0 6la ). The 
adiabatic functions \4>f dl ) are double-valued functions of 


= cos ^ I <j)f ,a ) + Sin I 0fj la ), (5) 

= - sin ^ | la ) + cos 10^ la ), (6) 
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the nuclear parameters (x, y) because encircling the Cl 
point corresponds to change in 9(x, y ) from 0 to 2tt which 
leads to a sign change in |</>f dl ). The unitary transforma¬ 
tion to the adiabatic representation brings i? 2 D to a form 


adi I T 2D + Tn 


rjadi _ 

^2D — 


T21 


Tl2 

? 2 D + T22 , 


W x 0 

0 w 2 , 


( 8 ) 


where 


vc 1 , 2 =i(y 1 i + y 2 2)=F^\/^^2) 7 74^, (9) 


are the adiabatic potentials with the minus (plus) sign 
for Wi (' W 2 ), and 


Tij — h~dij ■ V ^ Dij, (10) 

are kinetic energy terms containing NACs 

dy = <^ adi | V<Af), D tj = (#*1 V 2 0f) (11) 


B. Mixed quantum-classical methods 


One of the most straightforward routes to the MQC 
methods is to take the classical limit for nuclear DOF in 
the system wave-function, then for a two-state problem 
within the adiabatic representation an electronic wave- 
function is given 

2 

^ e (r; R, t) = '^2c j {t)<j>f l (r;'R(t)). (16) 

1=1 

The time-dependent coefficients are obtained from pro¬ 
jection of the electronic time-dependent Schrodingier 
equation (TDSE) 


2 2 

ihJ2[c^f + oof] = £ CfcW4(R)(/>f (17) 
1=1 fc=1 


onto the electronic adiabatic basis <^f (r; R(t)) 



/ Wi(R) — *fidi2 • R.\ /ci 
f M 2 i ■ R W 2 (R) J ^c 2 


with V = (d/dx,d/dy). Substituting definitions of the 
adiabatic states into Eq. (10), Tij can be express as 

A H 2 Dii H 2 (x 2 + y 2 ) 

2 8(7 _1 a; 2 + yy 2 ) 2 ’ 

H(L Z -L Z ) 

T d 4 ( 7 _1 a; 2 + 7 y 2 ) ’ * ^ 

where L z = —ih(xd/dy— yd/dx) is the angular momen¬ 
tum operator, and the overhead arrows indicate opera¬ 
tor’s directionality. 

In quantum dynamics within the adiabatic representa¬ 
tion, the GP can be introduced as a position-dependent 
phase factor e *®(*i3/)/ 2 f or single-valued nuclear basis func¬ 
tions, where 0(x,y) is provided by Eq. (7). 21 This phase 
factor would change signs of nuclear wave-functions on 
encircling the CI. However, it is more convenient to use 
this factor as a unitary transformation of the adiabatic 
Hamiltonian: = e~ lS ^ 2 H 2 ( ^e lS ^ 2 , then for simulat¬ 

ing GP effects we can use a normal single-valued nuclear 
basis with the H^ Hamiltonian. Also this approach al¬ 
lows one to compare Hf and H^d to see what operator 
terms are responsible for introducing GP effects. The 
comparison reveals that the phase factor alters kinetic 
energy terms 

GP = h(L z - L z ) h 2 {x 2 + y 2 ) 

4 ( 7 - 1 x 2 + 7 y 2 ) 4 ( 7 -1 £ 2 + 7 y 2 ) 2 ’ 

CP = h{L z - L z ) _ h 2 (x 2 +y 2 ) 

b 4(7 _1 a: 2 + 7y 2 ) 4(7 _1 x 2 + 7?/ 2 ) 2 ’ 


7 = 


X 0 U>{ 


( 12 ) 


(13) 


Here, the chain rule for NACs is used 

(f(R(t))|f(R(l))) = R.d J , (19) 

Thus, due to orthogonality of the adiabatic states the 
system population of each electronic state is given by 
\ci(t)\ 2 and the population transfer is regulated by the 
off-diagonal elements — ihdi 2 ■ R in Eq. (18). 

The nuclear equations of motion (EOM) for the EF 
method are derived from the total energy conservation 
condition 


dE 

dt 


d_ 

dt 



H^c 

- n adi C 


= 0, 


( 20 ) 


with = SijWi( R), which leads to Newton’s EOM 

for nuclei 


where 


R a =- 


Cf^ladic 

dR n 


c trH (e) d 

L adi’ U ° 


( 21 ) 


d n — 


0 

^! di lS 


{ft 


adi I ^4 > \ 



( 22 ) 


Thus, the classical particle moves under an averaged force 
produced by involved PESs. This is even more obvious 
if one reformulates the nuclear EOM in the diabatic rep¬ 
resentation 


R a = — 


d 

(Vn 

Vj 2 V 

dRa 

1^21 

V 22 J 


(23) 


This reformulation can be done either starting from the 
beginning using the diabatic representation 


and thus changes probabilities of electronic transitions. 


ip e (r,t) = Ci(t)(f>f la, (r) + c 2 (t)0 dla ( r ), (24) 
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or only by applying the adiabatic-to-diabatic transfor¬ 
mation in Eq. (21). This invariance with respect to the 
electronic state representation is one of the advantages 
of the EF method that is not shared by the SH method. 

In the SH case, nuclear EOM are also in the Newto¬ 
nian form but they evolve on a single electronic surface. 
There is some freedom in defining individual electronic 
surfaces on which nuclear dynamics takes place with the 
only constraint of the energy conservation when a sur¬ 
face switch (hop) takes place . 14 This freedom of choice 
in the electronic surface prompted some works where the 
DBOC had been added to the adiabatic states . 23,24 The 
rational can be given if we consider the full quantum nu¬ 
clear equation obtained by projecting the full TDSE onto 
the adiabatic electronic basis 


ih 



( F 2 D + Tn 

T12 \ 

y T21 

P 2 D +t 22 J 





(25) 

(26) 


nates centered at the Cl: 



/ h 2 d d 
\ P dp P dp 



12 


+ P 


( COS if 

sin^? 


sin ip 

— COS (p 


(28) 


where L z = —ihd/dip. It is easy to verify that Hm 
commutes with the vibronic angular momentum oper¬ 
ator J = L z 1 2 + 7 j (T y , where cr y is the Pauli matrix. 
Eigenfunctions of J are 


{ip |ro)f a = e imv 



and 


H m)f a = e im * 


(- sin f 

V cos 2 


(29) 


(30) 


Grouping all diagonal potential-like terms involves the 
second order NACs —H 2 Da /2 which are functions of R, 
therefore, they can be added to the potential energy sur¬ 
faces Wi to create modified surfaces 

Wj(R) = Wj(R) - yAi(R). (27) 

Considering that the DBOC has the h 2 prefactor, its ad¬ 
dition may seem insignificant. However, Du( R) diverges 
at the Cl point [Eq. (12), and Fig. 2], and hence, in the 
vicinity of the Cl the DBOC cannot be neglected based 
on its relatively small values far from the Cl. 

Note that there are no h terms present in the force 
definition in Eq. (21), this shows completely classical na¬ 
ture of the nuclear EOM. Also, even though we work in 
the adiabatic representation, the second order derivative 
couplings do not appear in the working equations because 
the nuclear kinetic energy has not been considered as a 
quantum operator. Besides the reason of inconsistency in 
powers of F, introducing the DBOC into the EF method 
would break the invariance of this method with respect 
to the electronic state representation. 


C. Mexican hat model 


where m are half-integer eigenvalues. Functions (ip |TO) dla 
are single-valued as eigenfunctions of a general quantum- 
mechanical operator without external parameter depen¬ 
dence should be. 

Let us now transform both operators to the adiabatic 
representation. For this model, the 6 angle [Eq. (7)] of 
the unitary transformation [Eqs. (5) and ( 6 )] becomes 
the polar angle ip. The unitary transformation of Hm 
and J leads to 


o-adi _ 1 ( & 9 9 , * , , 

Hm -2 d~p P d~P + P 1 2 


Tadi 


(L z 1 2 2 (Ty\ 


+ 2p2 +P<T Z , 

(31) 

l z i 2 . 

(32) 


The transformation of the J eigenfunctions [Eqs. (29) 
and (30)] gives (y>|m) adl = (e* mv ,0)^ and (ip [m)^' = 
(0, e lrnip y, which seem as regular eigenfunctions of L z 
apart from their half-integer values of m. Thus, the 
eigenfunctions of J adl are double-valued functions. 

Since the commutation relations are the same in all 
representations, and J adl commute and have a com¬ 
mon system of eigenfunctions, hence, the eigenfunctions 
of H^ 1 [Eq. (31)] can be sought as 


The adiabatic potentials IF,; in Eq. (9) of the 2D LVC 
Hamiltonian acquires a cylindrical symmetry for 7 = 1 
[Eq. (12)]. This symmetry facilitates comparison of non- 
adiabatic transfer elements for different methods. Thus, 
we will consider in details the 7 = 1 case, also known as 
the Mexican hat model. 

For our analysis, it is convenient to write H 2 d [Eq. (2)] 
with wi = uj 2 = xq = 1 and A = 0 in polar ( p , tp) coordi- 


f (P,¥?IXiA 
\(p,¥> Ix2>y 



(33) 


In Eq. (33) the double-valuedness of adiabatic nuclear 
wave-functions is isolated in (<p\m) = e mip . One can 
turn GP effects “on” and “off” by changing to: half¬ 
integer values correspond to inclusion of the GP, whereas 
integer values mean the GP is neglected. 
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To perform comparative analysis of quantum methods 
with and without GP with the same set of integer to an¬ 
gular functions, we apply a gauge transformation of Mead 
and Truhlar, 21 e*^/ 2 , to all half-integer angular functions 
|to). This transformation changes the Hamiltonian 
by modifying kinetic energy terms as in the case of the 
general 2D LVC model [Eqs. (15) and (14)]. For the 
Mexican hat model we can separate angular and radial 
components for all kinetic energy terms which allows for 
more detailed analysis. To estimate to dependence of 
quantum transition probabilities without breaking sym¬ 
metry with respect to left and right rotation around the 
Cl point we consider sums of t \2 elements projected onto 
|±m) | fi) states. For the Mexican hat Hamiltonian with¬ 
out explicit GP terms [Eq. (31)] transition amplitudes 
are 


ri 2 to = 


h 

+ 2 


= H\r 


</l| H ^2 l m ) I/ 2 ) 

</i| {-m\ ^ I to) 1 / 2 ) 


1 


(/i| ^2 I/ 2 ) 


Adding the GP modifies transition elements as 

H 


r 12 P (l m IJ = 


(A IH Lz 2 ^ 2 l m > IA) 


(A I {~m\ L \ |-m) I/2) 


2 p 2 


TO — o + TO + 


(A 


h 2 


IA) 


1 ¥ 

(Al 2^1/2) , 

h 2 \m\ 

</i| 2^2 I/2) 


V 

ItoI = 0 


(34) 

(35) 


(36) 

(37) 

(38) 


Note that the only difference from adding the GP correc¬ 
tion is a transition enhancement for the to = 0 compo¬ 
nent in the presence of the GP (Fig. 3). 

For the classical treatment of nuclear motion in MQC 
methods, we note that the classical nuclear angular mo¬ 
mentum is conserved because of the cylindrical symme¬ 
try. The transition amplitude is given by 


not change the outcome t ^^ C (| to |) = | t ^ < ^ C ( to )|/2 + 
l r i 2 b C ( — to )I/2 = T i 2 b C ( TO )- ^ one neglects the differ¬ 
ence between the classical 1/p 2 term and its quantum 
analogue, then for to = 0 and are the same 

(Fig. 3). This is a result of a continuous nature of the 
classical angular momentum that after integrating over 
the range |m| £ [0,1] introduces a GP-like enhancement. 



FIG. 3. Transition amplitude T 12 as a function of \m\ for dif¬ 
ferent methods: quantum without the GP Eq. (35) (blue), 
quantum with the GP Eq. (38) (black), continuous classi¬ 
cal Eq. (39) (red), binned classical Eq. (40) (magenta), p- 
dependent parts are neglected in all cases. 


However, the solely angular dependence consideration 
raises a question whether the continuous character of the 
classical angular momentum will cause overestimation of 
the transition probability for high |m|’s (Fig. 3, \m\ > 1). 
The answer is negative, and it becomes obvious if we 
account for the p-dependence of T 12 terms. Although p 
and L z are independent, intuitively, it is clear that for 
a general trajectory, due to the centrifugal force, high 
to’ s have low 1/p 2 factors in Eq. (40). Therefore, we can 
approximate (to + 1 /2)/p 2 ~ m/p 2 for large |m|’s. To 
confirm this intuitive consideration we performed both 
classical and quantum simulations for the Mexican hat 
model. In quantum simulations the initial state is given 
by a Gaussian wave packet placed on the upper cone 


t“ qc = *M 21 -R=^, (39) 

where L z = xP y — yP x is the classical angular momen¬ 
tum. For comparison with the quantum results we will 
perform quasi-classical binning by integrating continuous 
L z = hm values between discrete values of m and to + 1 


MQC 

r 12,b 


( to ) = h 



mdm 


H(m + 1/2) 

V 


(40) 


The same contribution will appear if we consider the 
[—to, —to—1 ] range, thus the averaging of two results does 


X2{x,y,0) 



exp 


(x - x 0 ) 2 


(7 


2 

X 



(41) 


with widths cr x = a y = \/2. The classical counter¬ 
part is initiated from a Gaussian distribution for posi¬ 
tions and momenta obtained via the Wigner transform 
of X2{x,y,0). Average transition amplitudes split to to 
components in Fig. 4 confirm the conjecture of the radial 
component (1/p 2 ) reduction with increase of the angu¬ 
lar component (to). Moreover, these effects compensate 
each other consistently through the m series so that both 
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TABLE I. Parameters of the 2D effective LVC Hamiltonian, 
Eq. (2) for the three studied systems. 


Wi 

UJ 2 a c 

A 

7 

7.743 x 10 -3 

Bis (methylene) adamantyl cation 
6.680 x 10 -3 31.05 8.092 x 10 -5 

0.000 

0.09 

9.557 x 10“ 3 

Butatriene cation 

3.3515 x 10~ 3 20.07 6.127 x 10 -4 

0.020 

0.67 

3.650 x 10“ 3 

Pyrazine 

4.186 x 10 -3 48.45 4.946 x 10~ 4 

0.028 

1.50 


methods plateaux at large |m|’s. Therefore, a continuous 
character of classical angular momentum helps to mimic 
the enhancement of the fully cylindrical m = 0 compo¬ 
nent and does not interfere with other angular compo¬ 
nents. 

2 - 1 - 1 - 1 - 1 - 1 - 1 - 


(!) 


Cl 

E 

ra 

c 

.9 

c 

2 

<13 

CO 

(13 

o 0.5 
> 

< 

0 t-- 1 - 1 - 1 - 1 - 1 - 

0 1 2 3 4 5 6 7 

M 

FIG. 4. Transition amplitudes, h 2 \m (/i| (2p) -2 I/ 2 ) | for 
quantum (sold red) and h(\m\ + l/2)(2p) -2 for MQC (dashed 
blue) simulations, averaged over the first a.u. of simulation 
time as a function of \m\. 


III. MOLECULAR CALCULATIONS 

Here we consider three molecular systems: the 
bis (methylene) adamantyl cation, (BMA ) 25 ’ 26 the buta- 
triene cation , 2 ' -32 and the pyrazine molecule . 33-35 These 
systems have been extensively studied before, and it was 
shown that they are well described by multi-dimensional 
LVC models. For our simulations, we have reduced N- 
dimensional LVC models to effective 2D LVC models us¬ 
ing collective nuclear DOF so that the 2D models repro¬ 
duce short-term dynamics of the original N-dimensional 
models [see the Appendix of Ref. 6 for details]. These 
three systems have LVC parameters (Table I) that are 
representative for various manifestations of GP effects . 6 

To provide comparative analysis of dynamical features 
appearing from DBOC and GP influences in addition to 


MQC simulations we provide quantum results obtained 
using three nuclear non-adiabatic Hamiltonians: 1) the 
full Hamiltonian [Eq. (2)] producing exact dynamics, 2) 
the “no GP” Hamiltonian [Eqs. ( 8 ), (12), and (13)], and 
3) the “no GP, no DBOC” Hamiltonian [Eqs. ( 8 ), (13), 
and Tu = 0]. Numerical procedures to propagate the 
TDSE with these Hamiltonians are detailed in Ref. 6 . 
For all methods we compare the adiabatic population of 
the excited electronic state (t) = (% 2 (t) |y 2 (Q), where 
X 2 (x, y, t) is a time-dependent nuclear wave-function that 
corresponds to the excited adiabatic electronic state ini¬ 
tiated as a Gaussian distribution in Eq. (41) with widths 
<j x = \j2/ui and a y = To calculate P^(t) in 

MQC simulations, |c 2 1 2 (^) are used in the EF approach, 
and the instantaneous ratios between the number of tra¬ 
jectories on the excited state to the total number of tra¬ 
jectories are taken in the SH approach. In both MQC 
methods, P^J (t) are averaged over 2000 trajectories with 
nuclear momenta and positions sampled from the Wigner 
transform of the corresponding quantum wave packet. To 
integrate MQC EOM, the 4th order Runge-Kutta inte¬ 
grator has been used with the time-step 0.05 fs for SH and 
0.001 fs for EF methods, respectively. For the SH method 
we used Tully’s fewest switches algorithm 14 with nuclear 
forces obtained from adiabatic PESs, W\ [Eq. (9)]. To 
illustrate the influence of the DBOC in SH dynamics, we 
introduce a modification, further referred as SH+DBOC, 
where nuclear forces are obtained from adiabatic PESs 
with the DBOC, W t [Eq. (27)]. 

As has been established in our previous study , 6 GP ef¬ 
fects manifest themselves quite differently in these molec¬ 
ular models: for BMA, GP effects are predominantly 
in compensation of a potential repulsion introduced by 
the DBOC, while for the butatriene cation and the 
pyrazine molecule the GP enhances non-adiabatic trans¬ 
fer of wave-packet’s m = 0 component. This difference 
stems from the anisotropy of the branching space, which 
is well characterized by the value of I 7 — 7 _1 |, the smaller 
this value is the closer the problem to the Mexican hat 
case is and more cylindrical all potential terms including 
the DBOC are. For BMA, 1 7 — 7 _1 | = 11.4 and this 
makes the DBOC a wide repulsive wall, while for the 
other systems, |y — 7 _1 | ~ 0.8 which is much closer to 
the Mexican hat limit (|y — 7 -1 | = 0). Thus, we will ini¬ 
tially analyze performance of the MQC methods for the 
BMA cation and then for the other two systems. 

a. BMA: DBOC in MQC. The exact population dy¬ 
namics in BMA demonstrates coherent oscillations that 
can be easily understood considering weak diabatic cou¬ 
plings in this system (Fig. 5). Thus, the exact dynam¬ 
ics almost solely undergoes on a single diabatic state 
that corresponds to the excited adiabatic state before the 
crossing and the ground adiabatic state after the cross¬ 
ing. Excited state populations in both MQC methods 
reproduce almost exactly those of the full QM dynamics. 
Smaller amplitudes of adiabatic population oscillations 
in SH than those in the exact dynamics is a manifesta¬ 
tion of SH overestimation of diabatic population trans- 
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0 5 10 15 20 25 

t, fs 

FIG. 5. Excited state adiabatic population () dynamics 
for the BMA cation obtained with different methods. 

fer. The origin of this overestimation and violation of 
the Marcus theory limit has been found in higher elec¬ 
tronic coherences within the SH approach. 36,37 The ex¬ 
planation of the overall success of the MQC methods for 
BMA is in the absence of both GP and DBOC terms 
in these methods. Therefore the MQC methods do not 
need the GP for cases where the main role of the GP is 
the DBOC compensation. On the other hand, compar¬ 
ing the SH+DBOC dynamics with the exact one shows 
that adding the DBOC to the adiabatic potentials can 
be very detrimental for MQC results (Fig. 5). The im¬ 
pact of uncompensated DBOC terms in MQC dynamics 
is even bigger than removing GP terms in quantum sim¬ 
ulations. This is a result of a repulsive potential of the 
DBOC that prevents classical nuclear dynamics in the SH 
method to approach a region of strong non-adiabatic cou¬ 
pling. Quantum wave packets for the “no GP” Hamilto¬ 
nian can increase the non-adiabatic transfer due to some 
tunnelling under the DBOC potential. Thus these results 
demonstrate that the DBOC should not be added to the 
MQC methods. 

Although the BMA branching space is significantly 
anisotropic, the transfer of the wave-packet m = 0 com¬ 
ponent is still somewhat suppressed. 6 Since the weight of 
the wave-packet m = 0 component near the Cl for BMA 
is quite substantial, 42%, absence of the GP enhancement 
of its transition in the “no DBOC, no GP” model results 
in deviation of the corresponding population dynamics 
from the exact one after 4 fs. Thus, even in anisotropic 
systems the GP induced non-adiabatic transfer enhance¬ 
ment and its imitation by MQC methods can be crucial 
for quantitative agreement with exact dynamics. 

b. and pyrazine: m = 0 enhancement in MQC. 

In the butatriene cation and the pyrazine molecule, the 
DBOC potential is relatively isotropic and compact. 
Therefore, it does not prevent a nuclear wave packet 
from accessing the vicinity of the Cl, and the DBOC 
presence does not change quantum dynamics significantly 
(see Figs. 6 and 7). Yet, GP effects are significant, be¬ 


cause the GP related term in accelerates the non- 
adiabatic transfer for the m = 0 component. Interest¬ 
ingly, this acceleration is well mimicked by the MQC 
methods due to the classical description of the angular 
momentum which leads to similar enhancement of the 
to = 0 component transfer. According to the angular 
decomposition of quantum wave packets at the moment 
of the closest proximity to the Cl, in both systems, the 
weight of the to = 0 component is close to 90%. There¬ 
fore, this enhancement is the main GP effect in these sys¬ 
tems. The excited state population dynamics in Figs. 6 
and 7 reveal that the MQC methods can reproduce the 
exact quantum dynamics and perform better than quan¬ 
tum methods that do not account for the GP. In both 
systems, the SH method performs slightly better com¬ 
pared to the EF approach. 



FIG. 6. Excited state adiabatic population (p£j?) dynamics 
of the butatriene cation obtained with different methods. 



FIG. 7. Excited state adiabatic population (F^di) dynamics 
of pyrazine obtained with different methods. 

















IV. CONCLUSIONS 

It has been recently found that pure quantum effects 
associated with the nontrivial geometric phase appearing 
in the nuclear and electronic adiabatic wave-functions for 
surface crossing problems can significantly affect popu¬ 
lation transfer dynamics. Although MQC methods ig¬ 
nore nuclear GP effects by substituting quantum nu¬ 
clear dynamics with its classical approximation, they are 
still very successful in simulating non-adiabatic dynamics 
through CIs. In this work we have unraveled the key ele¬ 
ments of this success: Both types of GP effects involved 
in the excited state dynamics, the DBOC compensation 
and the enhancement of the non-adiabatic transfer for the 
fully cylindrical component of a wave-packet, are mim¬ 
icked fortuitously in MQC methods using classical me¬ 
chanics. 

Interestingly, the DBOC term did not appear in orig¬ 
inal derivations of MQC schemes and have been added 
only later in an ad hoc manner. This work clearly demon¬ 
strated that such addition can be very detrimental for the 
quality of results and should be avoided. The mechanism 
for the cylindrical component enhancement in the MQC 
schemes has been elaborated on the Mexican hat model. 
The situation in some sense is opposite to the famous 
Planck quantization via discrete summation to describe 
the black-body radiation. In MQC transfer element, the 
purely quantum effect from the GP is recovered because 
a discrete summation over the angular eigenstates of the 
angular momentum operator is substituted by a classical 
continuous integration. Thus, nuclear GP effects make 
excited state non-adiabatic dynamics more classical by 
compensating some other quantum effects. 

There have been several proposals on using the 
Landau-Zener (LZ) formula 38,39 to model non-adiabatic 
dynamics through the conical intersection. 40-43 These 
proposals involve application of the LZ equation for prob¬ 
ability transfer with a subsequent averaging over individ¬ 
ual classical trajectories. Surprisingly, a question of in¬ 
fluence of the conical intersection topology on the result 
has not ever been raised. The current work can be used 
to rationalize an application of LZ-based approaches to 
the Cl problem even though in such methods topological 
geometric phase effects are not explicitly accounted. 
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